#########################################################################################
rm(list=ls(all=TRUE))
#install.packages("Amelia")
library(Amelia)
#install.packages("readstata13")
library(readstata13)
#install.packages("Hmisc")
library(Hmisc)
#install.packages("maptools")
library(maptools)

# set working directory
#setwd()

mydata <- read.dta13("IndustryLevelFDI.dta")

#########################################################################################
###################### Generate 10 imputed datasets for each industry ###################

for (i in 1:18){
mydata1 <- mydata[,c(1:19,19+i)]
head(mydata1)

set.seed(231011); imputation <- amelia(mydata1,m=10,ts="year",cs="country",polytime=3)
write.amelia(obj=imputation,file.stem=paste("Industry",i,"",sep="_"),row.names=F)
}

#########################################################################################
#########################################################################################
#Execute  savecoefs.do in Stata to retrieve coefficiet estimates for each imputed dataset
#########################################################################################
#########################################################################################


#########################################################################################
########### combine coefficients estimates from multiple imputed datasets ###############
rm(list=ls(all=TRUE))
mydata <- read.dta13("coefs.dta")
coef.matrix<-matrix(NA,15,36)

for (i in 1:18){
mydata1 <- subset(mydata,industry==i)
coefs<- mydata1[,c("var","coef","imputation")]
coefs <- reshape(coefs,timevar = "var",idvar="imputation",direction = "wide")
coefs <- coefs[,-1]

se <- mydata1[,c("var","stderr","imputation")]
se <- reshape(se,timevar = "var",idvar="imputation",direction = "wide")
se <- se[,-1]
result <- mi.meld(q = coefs, se = se)
coef.matrix[,c(2*i-1,2*i)]<-matrix(cbind(result$q.mi, result$se.mi),15,2 )
}
rownames(coef.matrix)<-unique(mydata[,1])

write.csv(coef.matrix,file = "results_allindustries.csv",row.names = T)


#########################################################################################
######################## Function: Coef Plot (lines 60-81) ##############################
coefplot<-function(x,coefs,SEs,adj.x,labelname){
  par(mfcol = c(1, 1), mar=c(3.5,3.5,0.5,1), cex = .6, cex.axis = .8)
  ci95.u<-coefs+1.96*SEs
  ci95.l<-coefs-1.96*SEs
  ci90.u<-coefs+1.65*SEs
  ci90.l<-coefs-1.65*SEs
  
  plot(x,coefs,pch=16, 
       ylim=c(min(ci95.l,na.rm=T)-0.02,max(ci95.u,na.rm=T)+0.02),
       xlim=c(min(x,na.rm=T)+adj.x,max(x,na.rm=T)),
       xlab="Capital Expenditures/Output (log)",ylab="Estimated Coefficients",
       mgp = c(1.8,.6,0),cex.lab=.8, cex.axis=.6,cex=.7)
  
  abline(lm(coefs~x),lwd=0.8,lty=2)
  abline(h=0,lty=2, col="gray",lwd=0.8)
  
  segments(x,ci95.l, x,ci95.u,lwd=0.6)
  segments(x,ci90.l, x,ci90.u,lwd=1.2)
  
  pointLabel(x,coefs+0.01,paste(labelname), cex=.6, allowSmallOverlap=F,
             offset=-0.05, pos=3)
}
########################################################################################
################################# Coef Plot: Figure 8 ##################################

dat.coef<-read.csv("results_allindustries.csv",header=T)
dat.coef1<-dat.coef[,-1]
coef.index<-seq(1,35,2)
se.index<-seq(2,36,2)
coef.per<-dat.coef1[1,coef.index]
se.per<-dat.coef1[1,se.index]

dat<-read.dta13("FixedAssetIntensity_US.dta")

dat$coefs<-t(coef.per)
dat$se<-t(se.per)

png(file="Personal-IndustryFDI-Imputed.png", width=4.5,height=3,unit="in",res=600)
coefplot(log(dat$FixAssetIntensity),dat$coefs,dat$se,adj.x=-0.09,labelname=dat$industry)
dev.off()


#########################################################################################
#########################################################################################
#Execute savecoefs.do in Stata to retrieve coefficient estimates for each industry(non-imputed)
#########################################################################################
#########################################################################################

#########################################################################################
################################# Coef Plot: Figure G-1 ###############################
dat.coef2<-read.dta13("coefs_nonimputed.dta")
dat.coef2$indcode2<-18:1
dat.coef2a<-dat.coef2[order(dat.coef2$indcode2),]

png(file="Personal-IndustryFDI-Nonimputed.png", width=4.5,height=3,unit="in",res=600)
coefplot(log(dat$FixAssetIntensity),dat.coef2a$coef,dat.coef2a$stderr,adj.x=-0.1,
         labelname=dat$industry)
dev.off()

#########################################################################################
################################ Table G-1 ############################################
coefs.out<-matrix(coef.matrix[1,],2,18)

coefs.m<-round(cbind( dat.coef2a$coef,dat.coef2a$stderr,
                      dat.coef2a$N,coefs.out[1,],coefs.out[2,],1782),digits=3)
coefs.m<-cbind(coefs.m, dat[,2],round(dat[,3],digits=3))

colnames(coefs.m)<-c("Coef1","SE1","N1","Coef2","SE2","N2",
                     "Industry","FixAssetIntensity")

write.csv(coefs.m[order(coefs.m[,8]),],file = "coefs_personalism.csv",row.names = T)
